In [1]:
cd C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents
In [3]:
import geopandas as gpd
import matplotlib.pyplot as plt
import pickle
import pandas as pd
import numpy as np
import calendar
In [2]:
import matplotlib.cm as cm
import matplotlib.colors as cls
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
In [3]:
#latex
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'custom'
matplotlib.rcParams['mathtext.rm'] = 'Times New Roman'
matplotlib.rcParams['mathtext.it'] = 'Times New Roman:italic'
matplotlib.rcParams['mathtext.bf'] = 'Times New Roman:bold'
In [4]:
from scipy.stats import boxcox, probplot
In [4]:
bios = gpd.read_file('biomes_hemispheres')
In [5]:
bios.columns
Out[5]:
In [7]:
area = bios['SUM_Area']*1e-6
In [10]:
data_columns = bios.columns[3:3+17*12]
In [11]:
inter_cols = {}
for yr in range(17):
inter_cols[yr] = data_columns[12*yr:12*yr +12]
In [12]:
intra_cols = {}
for mo in range(12):
intra_cols[mo] = [data_columns[mo + yr *12 ] for yr in range(17)]
for yr in range(17):
inter = bios[ inter_cols[yr]].sum(axis=1)/area
label = 'inter_' + str(yr +2001)
bios= bios.assign(label=inter)
bios = bios.rename(columns={'label': label})
inter_columns = bios.columns[-17:]
inter_max = bios[inter_columns].max().max()
for mo in range(12):
intra = bios[ intra_cols[mo]].sum(axis=1)/(17*area)
label = 'intra_{:02}'.format(mo)
bios = bios.assign(label=intra)
bios = bios.rename(columns={'label': label})
intra_columns = bios.columns[-12:]
In [18]:
intra_max = bios[intra_columns].max().max()
In [76]:
def draw_map(df,column,title=None,vmin=0,vmax=1,name='name',dpi=100,cmap='OrRd'):
fig, ax = plt.subplots(figsize=(12, 8))
ax.set_aspect('equal')
ax = df.plot(ax =ax, column=column, cmap=cmap,vmin=vmin, vmax=vmax)
ax.set_axis_off()
norm = cls.Normalize(vmin, vmax)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable._A = []
cbaxes = inset_axes(ax, width="80%", height="1%", loc=8)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig(name, dpi=dpi, bbox_inches='tight')
In [ ]:
''' backup
def draw_map(df,column,title=None,vmin=0,vmax=1,name='name',dpi=100):
fig, ax = plt.subplots(figsize=(12, 8))
ax.set_aspect('equal')
ax = df.plot(ax =ax, column=column, cmap='OrRd',vmin=vmin, vmax=vmax)
ax.set_axis_off()
norm = cls.Normalize(vmin, vmax)
cmmapable = cm.ScalarMappable(norm, 'OrRd')
cmmapable._A = []
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig(name, dpi=dpi, bbox_inches='tight')
'''
Inter-anual¶
In [78]:
lambdas = []
for yr in range(17):
column= 'inter_{}'.format(2001 + yr)
test_column, lbd = boxcox(bios[column] + 1e-25)
lambdas.append(lbd)
lmbda= sum(lambdas)/len(lambdas)
In [79]:
lambdas
Out[79]:
In [22]:
probplot(bios.inter_2017, plot=plt.gca());
In [23]:
probplot(test_column, plot=plt.gca());
In [80]:
maxs = []
mins = []
for yr in range(17):
column= 'inter_{}'.format(2001 + yr)
test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
maxs.append(test_column.max())
mins.append(test_column.min())
box_min = min(mins)
box_delta =max(maxs) - box_min
In [81]:
for yr in range(17):
column = 'inter_{}'.format(yr +2001)
test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
bios = bios.assign(test_column =test_column)
yr_name = str(yr +2001)
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome in each hemisphere ' +yr_name
name = 'inter_biomes_nsparts_index_{}'.format(yr +2001)
draw_map(bios,'test_column',title=title,name=name)
Intra-anual¶
In [82]:
lambdas = []
for mo in range(12):
column= 'intra_{:02}'.format(mo)
test_column, lbd = boxcox(bios[column] + 1e-25)
lambdas.append(lbd)
lmbda= sum(lambdas)/len(lambdas)
In [84]:
lambdas
Out[84]:
In [83]:
probplot(bios.intra_11, plot=plt.gca());
In [29]:
probplot(test_column, plot=plt.gca());
In [85]:
maxs = []
mins = []
for mo in range(12):
column= column= 'intra_{:02}'.format(mo)
test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
maxs.append(test_column.max())
mins.append(test_column.min())
box_min = min(mins)
box_delta =max(maxs) - box_min
In [86]:
for mo in range(12):
column = 'intra_{:02}'.format(mo)
test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
bios= bios.assign(test_column =test_column)
mo_name = calendar.month_name[mo + 1]
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome in each hemisphere ' +mo_name + ' (2001-2017)'
name = 'intra_biomes_nsparts_index_{:02}'.format(mo)
draw_map(bios,'test_column',title=title,name=name)
Total¶
In [32]:
total = bios[data_columns].sum(axis=1)/17
total = total/area
bios= bios.assign(total=total)
In [77]:
column = 'total'
test_column, lbd = boxcox(bios[column] + 1e-25)
test_column = (test_column - test_column.min())/(test_column.max() - test_column.min())
bios = bios.assign(test_column =test_column)
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome in each hemisphere (2001-2017)'
name = 'total_biomes_nsparts_index'
draw_map(bios,'test_column',title=title,name=name,dpi=600)
Anomalias¶
Inter-anual¶
In [38]:
inter_medias = bios[inter_columns].mean(axis=1)
In [41]:
matriz = np.array([[4,6,8],[6,8,10]])
In [43]:
matriz - np.array([[0],[2]])
Out[43]:
In [53]:
matriz.shape
Out[53]:
In [54]:
np.array([[0],[2]]).shape
Out[54]:
In [58]:
anomalias = np.array(bios[inter_columns]) - inter_medias.values.reshape(27,1)
In [45]:
bios[inter_columns]
Out[45]:
In [61]:
anomalias[0,:]
Out[61]:
In [56]:
inter_medias.values.reshape(27,1)
Out[56]:
In [64]:
anomalias.max()
Out[64]:
In [65]:
for yr in range(17):
test_column = anomalias[:,yr]
bios = bios.assign(test_column =test_column)
yr_name = str(yr +2001)
#"Burnt area pixels (#/km2) by biome in 2001
title = r'Anomalies burned pixels (#/$\mathrm{km}^{2}$) by biome N-S parts in ' +yr_name
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'inter_anomalies_biomes_nsparts_index_{}'.format(yr +2001)
draw_map(bios,'test_column',title=title,name=name,cmap='seismic',vmin=anomalias.min(),vmax=anomalias.max())
Intra-anual¶
In [87]:
intra_medias = bios[intra_columns].mean(axis=1)
In [88]:
anomalias = np.array(bios[intra_columns]) - intra_medias.values.reshape(27,1)
In [89]:
bigger = abs(anomalias.max()) if abs(anomalias.max()) >= abs(anomalias.min()) else abs(anomalias.min())
vmax = bigger; vmin = -1 * bigger
In [91]:
for mo in range(12):
test_column = anomalias[:,mo]
bios= bios.assign(test_column =test_column)
mo_name = calendar.month_name[mo + 1]
#title = r'Anomalies burned area pixels (#/$\mathrm{km}^{2}$) by biome in each hemisphere ' +mo_name + ' (2001-2017)'
title = r'Anomalies (#/$\mathrm{km}^{2}$) by biome ' +mo_name + ' (2001-2017)'
name = 'intra_anomalies_biomes_nsparts_index_{:02}'.format(mo)
draw_map(bios,'test_column',title=title,name=name,cmap='seismic',vmin=vmin,vmax=vmax,dpi=600)
In [93]:
anomalias.shape
Out[93]:
Total¶
In [73]:
anomalia = bios.total -bios.total.mean()
In [74]:
anomalia
Out[74]:
In [75]:
test_column = anomalia
bios= bios.assign(test_column =test_column)
title = r'Anomalies burned area pixels (#/$\mathrm{km}^{2}$) by biome N-S parts in (2001-2017)'
name = 'total_anomalies_biomes_nsparts'
draw_map(bios,'test_column',title=title,name=name,cmap='seismic',vmin=anomalia.min(),vmax=anomalia.max())
In [94]:
def test_map(df,column,title=None,vmin=0,vmax=1,name='name',dpi=100,cmap='OrRd'):
fig, ax = plt.subplots(figsize=(12, 8))
ax.set_aspect('equal')
ax = df.plot(ax =ax, column=column, cmap=cmap,vmin=vmin, vmax=vmax)
ax.set_axis_off()
norm = cls.Normalize(vmin, vmax)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable._A = []
cbaxes = inset_axes(ax, width="80%", height="1%", loc=8)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
fig.tight_layout()
plt.show()
fig.savefig(name, dpi=dpi, bbox_inches='tight')
In [95]:
test_map(bios,'test_column',title='teste',name='teste_tight',cmap='seismic',vmin=vmin,vmax=vmax,dpi=600)
Análise sem mapas¶
In [6]:
bios.columns
Out[6]:
In [7]:
len(bios.columns)
Out[7]:
In [10]:
relevant_columns = bios.columns[1:-3]
In [11]:
relevant_columns
Out[11]:
In [12]:
bios = bios [relevant_columns]
In [13]:
type( bios)
Out[13]:
In [64]:
bios
Out[64]:
In [15]:
from pandas import ExcelWriter
import xlsxwriter
In [16]:
#writer = ExcelWriter('results_newstyle2.xlsx',engine='xlsxwriter')
#df_finall.to_excel(writer,sheet_name='Burned Pixels',engine='xlsxwriter')
#writer.save()
In [21]:
bios['SUM_Area'] = bios['SUM_Area'].apply (lambda x:x/1e6)
In [25]:
bio_names = pd.read_excel('biomas_nomes.xlsx')
In [27]:
bio_names.set_index('BIOME_NUM')
Out[27]:
In [35]:
type(bio_names.loc[8].values[1])
Out[35]:
In [36]:
def bioname(bionum):
return bio_names.loc[bionum].values[1]
In [37]:
bionames = bios['BIOME_NUM'].apply(lambda x:bioname(x))
In [39]:
bios = bios.assign(biome_names = bionames)
In [41]:
bios.columns
Out[41]:
In [45]:
new_columns =[ bios.columns[0]] +[ bios.columns[-1]] + list(bios.columns[1:-1] )
In [46]:
new_columns
Out[46]:
In [47]:
bios = bios[new_columns]
In [51]:
data_columns = bios.columns[3:-1]
In [52]:
data_columns
Out[52]:
In [49]:
writer = ExcelWriter('biome_analysis.xlsx',engine='xlsxwriter')
bios.to_excel(writer,sheet_name='original',engine='xlsxwriter', index= False)
In [53]:
inter_cols = {}
for yr in range(17):
inter_cols[yr] = data_columns[12*yr:12*yr +12]
In [ ]:
intra_cols = {}
for mo in range(12):
intra_cols[mo] = [data_columns[mo + yr *12 ] for yr in range(17)]
In [58]:
area = bios['SUM_Area']
In [59]:
for yr in range(17):
inter = bios[ inter_cols[yr]].sum(axis=1)/area
label = 'inter_' + str(yr +2001)
bios= bios.assign(label=inter)
bios = bios.rename(columns={'label': label})
In [61]:
inter_columns = bios.columns[-17:]
In [67]:
inter_columns
Out[67]:
In [63]:
for mo in range(12):
intra = bios[ intra_cols[mo]].sum(axis=1)/area
label = 'intra_{:02}'.format(mo)
bios = bios.assign(label=intra)
bios = bios.rename(columns={'label': label})
In [65]:
intra_columns = bios.columns[-12:]
In [66]:
intra_columns
Out[66]:
In [57]:
intra_cols
Out[57]:
In [71]:
total = bios[data_columns].sum(axis=1)
In [73]:
total_area = bios[data_columns].sum(axis=1)/area
In [75]:
colunas_hemis = list(bios.columns[:3]) + ['SUM_Area']
In [76]:
colunas_hemis
Out[76]:
In [77]:
table_hemis = bios[colunas_hemis]
In [79]:
table_hemis = table_hemis.assign(total = total)
In [80]:
table_hemis = table_hemis.assign(total_area = total_area)
In [82]:
table_hemis['biome_names'][11]
Out[82]:
In [83]:
table_hemis.to_excel(writer,sheet_name='totais',engine='xlsxwriter', index= False)
In [84]:
writer.save()
In [85]:
table_hemis
Out[85]:
In [86]:
table_hemis.columns
Out[86]:
In [92]:
dict_agg = {'is_north': 'first', 'biome_names': 'first', 'SUM_Area': 'sum', 'total':'sum'}
In [96]:
table_hemis = table_hemis[['BIOME_NUM', 'biome_names', 'is_north', 'SUM_Area', 'total']].groupby('BIOME_NUM').agg(dict_agg).drop('is_north', axis=1)
In [98]:
total_area= table_hemis['total']/ table_hemis['SUM_Area']
In [99]:
table_hemis = table_hemis.assign(total_area = total_area)
In [108]:
table_hemis
Out[108]:
In [109]:
writer = ExcelWriter('biome_analysis2.xlsx',engine='xlsxwriter')
table_hemis.to_excel(writer,sheet_name='totais2',engine='xlsxwriter')
In [110]:
writer.save()
In [112]:
from openpyxl import load_workbook
In [135]:
wb = load_workbook('biome_analysis.xlsx')
In [136]:
writer = pd.ExcelWriter('biome_analysis.xlsx',engine='openpyxl')
In [137]:
writer.book = wb
In [138]:
writer.sheets = dict((ws.title, ws) for ws in wb.worksheets)
In [139]:
table_hemis.to_excel(writer, sheet_name='total_2')
In [140]:
writer.save()
In [142]:
#código acima, supostamente abria ficheiro excel já existente, e acrescentava uma nova folha. na verdade, acrescenteu uma
#nova folha mas não ao ficheiro já existente
In [143]:
bios
Out[143]:
In [148]:
bios[inter_columns]['inter_2001'].plot(kind='bar')
Out[148]:
In [150]:
bios.columns[:3]
Out[150]:
In [163]:
labels = ['{} {}'.format( int(bios['BIOME_NUM'].iloc[i]), bios['is_north'].iloc[i].upper()) for i in bios.index ]
In [164]:
labels
Out[164]:
In [175]:
ax = bios[inter_columns].plot(kind='bar')
ax.set_xticklabels( labels);plt.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
#fig = plt.gcf()
plt.figure(figsize=(16,24))
plt.show()
In [159]:
ax
Out[159]:
In [160]:
ax.set_xlabel =labels;plt.show()
inter_anual¶
In [184]:
fig, ax = plt.subplots(figsize=(24,16))
bios[inter_columns].plot(kind='bar',ax= ax)
ax.set_xticklabels( labels, rotation =None,fontsize=15)
plt.legend(fontsize=15)
plt.show()
In [ ]:
fig, ax = plt.subplots(figsize=(24,16))
bios[inter_columns].plot(kind='bar',ax= ax)
ax.set_xticklabels( labels, rotation =None,fontsize=15)
plt.legend(fontsize=15)
plt.show()
In [211]:
#fig, axes =plt.subplots(3,6,figsize =(24,16))
bios[inter_columns].plot(kind='bar',subplots=True, grid=True,title="Variação interanual", layout=(6, 3), sharex=True, sharey=True, legend = False, )
[ax.set_xticklabels(labels) for ax in plt.gcf().axes]
plt.tight_layout()
plt.show()
In [241]:
'''
fig, axes =plt.subplots(6,3,figsize =(24,16))
bios[inter_columns].plot(axes =axes ,kind='bar',title="Variação interanual", legend = False, )
[ax.set_xticklabels(labels) for ax in plt.gcf().axes]
plt.tight_layout()
plt.show()
'''
Out[241]:
In [240]:
'''
fig, axes = plt.subplots(nrows=6, ncols=3)
for i, c in enumerate(bios[inter_columns].columns):
bios[inter_columns][c].plot(kind='bar', ax=axes[i])
plt.tight_layout()
plt.show()
'''
Out[240]:
In [256]:
for i, c in enumerate(bios[inter_columns].columns):
bios[inter_columns][c].plot(kind='bar', title = c[-4:])
plt.gca().set_xticklabels( labels)
plt.show()
In [223]:
# não fucinona por causa das areas
#bios[list(inter_columns) + ['BIOME_NUM']].groupby('BIOME_NUM').sum()
Intra-anual¶
In [224]:
fig, ax = plt.subplots(figsize=(24,16))
bios[intra_columns].plot(kind='bar',ax= ax)
ax.set_xticklabels( labels, rotation =None,fontsize=15)
plt.legend(fontsize=15)
plt.show()
In [ ]:
month_names
In [230]:
for i, c in enumerate(bios[intra_columns].columns):
bios[intra_columns][c].plot(kind='bar', title = calendar.month_name[int(c[-2:]) +1][:3])
plt.gca().set_xticklabels( labels)
plt.show()
In [229]:
calendar.month_name[int('intra_00'[-2:]) +1][:3]
Out[229]:
In [232]:
bios.columns[:-17-12]
Out[232]:
In [233]:
biomes = bios[bios.columns[:-17-12]]
In [238]:
biomes = biomes.groupby('BIOME_NUM').sum()
In [239]:
biomes.columns
Out[239]:
In [242]:
inter_cols_t = {}
for yr in range(17):
inter_cols_t[yr] = biomes.columns[12*yr:12*yr +12]
In [244]:
area_t = biomes['SUM_Area']
In [246]:
for yr in range(17):
inter = biomes[ inter_cols[yr]].sum(axis=1)/area_t
label = 'inter_' + str(yr +2001)
biomes= biomes.assign(label=inter)
biomes = biomes.rename(columns={'label': label})
In [249]:
inter_columns_t = biomes.columns[-17:]
In [250]:
inter_columns_t
Out[250]:
In [254]:
fig, ax = plt.subplots(figsize=(24,16))
biomes[inter_columns_t].plot(kind='bar',ax= ax)
ax.set_xticklabels( range(15), rotation =None,fontsize=15)
plt.legend(fontsize=15)
plt.show()
In [255]:
for i, c in enumerate(biomes[inter_columns_t].columns):
biomes[inter_columns_t][c].plot(kind='bar', title = c[-4:])
#plt.gca().set_xticklabels( labels)
plt.show()
In [263]:
biomes[inter_columns_t].iloc[0].plot(kind='bar', title =table_hemis.iloc[0]['biome_names']); plt.gca().set_xticklabels( range(2001,2018))
Out[263]:
In [261]:
table_hemis.iloc[0]['biome_names']
Out[261]:
In [264]:
for i in range(15):
biomes[inter_columns_t].iloc[i].plot(kind='bar', title =table_hemis.iloc[i]['biome_names'])
plt.gca().set_xticklabels( range(2001,2018))
plt.show()
In [265]:
intra_cols_t = {}
for mo in range(12):
intra_cols_t[mo] = [biomes.columns[mo + yr *12 ] for yr in range(17)]
In [266]:
intra_cols_t
Out[266]:
In [ ]:
for mo in range(12):
intra = bios[ intra_cols[mo]].sum(axis=1)/area
label = 'intra_{:02}'.format(mo)
bios = bios.assign(label=intra)
bios = bios.rename(columns={'label': label})
In [267]:
#calendar.month_name[int(c[-2:]) +1][:3]
In [268]:
months_names = [calendar.month_name[i +1][:3] for i in range(12)]
In [282]:
for i in range(bios.shape[0]):
title ='{} {}'.format(table_hemis.iloc[int(bios.iloc[i]['BIOME_NUM'])]['biome_names'],bios.iloc[i]['is_north'].upper())
bios[intra_columns].iloc[i].plot(kind='bar', title=title)
plt.gca().set_xticklabels( months_names)
plt.show()
In [272]:
bios.shape[0]
Out[272]:
In [276]:
i= 3
In [277]:
'{} {}'.format(table_hemis.iloc[bios.iloc[i]['BIOME_NUM']]['biome_names'],bios.iloc[i]['is_north'].upper())
In [278]:
table_hemis.iloc[bios.iloc[i]['BIOME_NUM']]['biome_names']
In [280]:
table_hemis.iloc[bios.iloc[i]['BIOME_NUM']]